Work fluctuations in a time-dependent harmonic potential: 
rigorous results and beyond the overdamped limit 
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We investigate the stochastic motion of a Brownian particle in the harmonic potential with a 
time-dependent force constant. It may describe the motion of a colloidal particle in an optical trap 
where the potential well is formed by a time-dependent field. We use the path integral formalism 
to solve the Langevin equation and the associated Fokker-Planck (Kramers) equation. Rigorous 
relations are derived to generate the probability density function for the time-dependent nonequi- 
librium work production beyond the overdamped limit. We find that the work distribution exhibits 
an exponential tail with a power-law prefactor, accompanied by an interesting oscillatory feature 
(multiple pseudo locking-unlocking transitions) due to the inertial effect. Some exactly solvable 
cases are also discussed in the overdamped limit. 

PACS numbers: 05.70.Ln, 05.10.Gg, 05.40.-a 



I. INTRODUCTION 



Nonequilibrium (NEQ) fluctuation has been an impor- 
tant issue in the field of statistical mechanics for the 
last two decades since the first discovery of the fluctu- 
ation theorem for the entropy production [iHl]. Fluctu- 
ation theorems (FT's) [4l-[ll| and the Jarzynski equality 
(JE) [13, [l^l are the central theoretical relations govern- 
ing NEQ fluctuation phenomena, widely valid for many 
NEQ systems, deterministic or stochastic, thermostat- 
ted with heat baths. It basically deals with a large 
fluctuation around the average measurement in theory 
and experiment with considerable contribution from rare 
events. Such phenomena become dominant for a system 
with small degrees of freedom. There have been exten- 
sive studies of small experimental systems such as a mi- 

lj|, a single 



croscopic bead dragging in a viscous fluid ^-^j, _ 

molecule of the RNA under mechanical stretch [isl . fig , 
an oscillatin g: be ad under the translating center of the 
optical trap [131 > the circuit of an electric dipole in elec- 
tric potential bias (isj , an ultra light metallic wire under 
torsion [l^, etc. 

External bias was considered as a typical underlying 
mechanism for NEQ systems such as the nano circuit de- 
vice with potential bias [H, [lO] , the harmonic oscillator 
under constant torque applied [l^, [2l|, [^^ , one dimen- 
sional lattice gas in contact at boundaries with different 
heat or particle baths [l^, [l^ . A nonconservativc force 
was also recognized as a source for the entropy produc- 
tion such as in a nano heat engine in contact with mul- 
tiple reservoirs for a circulating current in high dimen- 
sional systems j25l - l27j . Non-Markovian nature caused 
by memory effect or colored noise is another source for 
NEQ [2^, [2^. In these examples, the system reaches a 
NEQ steady state (NESS) after a transient period, where 
a persistent nonzero current, directed or circulated, gen- 
erates the incessant work production. The probability 



distribution function (PDF) of the work production ex- 
hibits an exponential decay with a power-law prefactor in 
the rare-event region [23] , along with interesting unusual 
features such as initial condition dependency of the large 
deviation function [30l - [33j and multiple dynamic transi- 
tions in reaching the NESS [HHl]. 

On the other hand, a time-dependent perturbation 
on external parameters such as electric field, magnetic 
field, volume, force constant, etc. generates a genuine 
time-dependent NEQ state where the system never main- 
tains the NESS. For example, the stochastic motion of a 
Brownian particle was studied in the harmonic potential 
moving with a constant velocity {sliding parabola poten- 
tial) [35l-t40l|. and also in the harmonic potential with 
a time- dependent force constant {breathing parabola po- 
tential) [4l| - [44| . In these cases, the work PDF also shows 
an exponential decay with a power-law prefactor in the 
rare-event region, along with a time-dependent charac- 
teristic value for the work production determining the 
exponential decay shape. Most of previous studies con- 
sidered the overdamped limit, partly because the exper- 
imental situation for a colloidal particle in a harmonic 
trap can be well approximated in the overdamped limit 
and also partly because the analytical treatment is much 
easier [35l - [44| . 

In this paper, we generalize these results beyond the 
overdamped limit {underdamped case) for a Brownian 
particle in a breathing parabola potential with the mo- 
mentum variable kept intact. We focus on the inertial 
effect on the time-dependent characteristic value for the 
work production. Our model may also serve as a soft-wall 
version of the box expansion or compression with a sin- 
gle Brownian particle inside, in contact with a thermal 
reservoir [4^. The experimental setup is also feasible: 
In a molecular tweezer or an optical trap experiment, 
the potential well can be approximated by the harmonic 
potential. The shape of the potential well, so the force 
constant of the harmonic potential is set to vary with a 
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time-dependent external field. 

The stochastic motion is described by the Langevin 
equation and the corresponding Fokkcr-Planck 
(Kramers) equation. We use the path integral for- 
malism to derive rigorous relations from which the 
time-dependent work PDF can be easily calculated 
with machine accuracy and also its cumulants at all 
temperatures with any time-dependent force constant. 
We find an interesting oscillatory feature of the work 
PDF shape, solely due to the inertial effect (absent 
in the overdampcd limit), which resembles multiple 
dynamic transitions found in the linear diffusion sys- 
tem [23, [sj], but shows smooth crossovers rather than 
sharp transitions. Thus we call this crossover as a 
pseudo locking-unlocking transition. The existence of 
multiple pseudo dynamic transitions may be related to 
the existence of the phase-space circulating current in 
the underdampcd case, but full intuitive understanding 
calls for further investigation in future. 

In section ini we introduce the breathing harmonic po- 
tential function and discuss the FT's. In section Hill we 
derive the equations for the PDF and the cumulants for 
the work production using the path integral formalism. 
Our formalism is tested for the systems in the sudden 
change limit. In section ITVl we present the analysis for 
the work PDF and find the exponential tail with a power- 
law prcfactor. In section |Vl we study the overdampcd 
limit for exactly solvable cases. In section I VII we dis- 
cuss the main results of our work and the perspective to 
future works. 



II. TIME-DEPENDENT HARMONIC 
POTENTIAL 



We consider the Brownian motion of a particle in one 
dimension under the breathing harmonic potential with 
a time-dependent force constant k = k(t) and in contact 
with a heat bath. The equations of motion are given by 



X = p/m 

p = —jp/m. — kx + ^ 



(1) 



where 7 is a damping coefhcient and ^ is white noise 
with zero mean satisfying {^{t)£,{t')) = 2d S{t — t'). The 
diffusion coefficient d is chosen to satisfy the Einstein 
relation, d ~ /3~^7; which guarantees the equilibrium 
(EQ) Boltzmann distribution at inverse temperature j3 
in the steady state, if k is constant in time. 
The equations of motion can be rewritten as 



Fq + C, 



(2) 



where q = {x,p)'^ and ^ = (0,^)-^. Here, the superscript 
T denotes the transpose of a vector or a matrix. The 
force matrix F is given by 



The energy of a particle is given by E{q- k) = + 
which is written as E = ^q-^ • H • q with a Hamiltonian 
matrix 



k 
1/m 



(4) 



Let P(q, t) be the probability density function for find- 
ing a particle at state q and at time t. Then it satisfies 
the Fokker-Planck equation, specially called the Kramers 
equation. 



dP{q,t) 
Ft 



= V-(F.q + D. V)P(q,0 , (5) 



where V = {dx, dpf" and the diffusion matrix is given by 



D 



e 
d 



(6) 



where a small positive parameter e is introduced for con- 
venience, making possible the inversion of the diffusion 
matrix D during a formal manipulation in the path in- 
tegral formulation. In the end, we take the e limit 
to recover the delta function constraint 5{x — p/m) for 
position and momentum. 

With e, position and momentum can be treated on the 
same footing, which gives us the formal advantage over 
the usual path integral with the (5- function constraint. 
This approach works well. For instance, one can repro- 
duce the expected results for the EQ PDF when k is 
a time-independent constant [1^. In this case, the EQ 
Boltzmann distribution 
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Z{k) 



(7) 



becomes the stationary solution of the Kramers equation 
in the limit e 0. The partition function is given by 
Z{k) / dq e-'3^(i''=' (47r2m/(/32/!c))i/2^ so that the 
free energy is given by J^{k) = — ^ \n{ATT^m/ {/3^k)). 

When the force constant k varies in time, the system 
is driven into a NEQ state. It belongs to the Jarzynski's 
criterion for NEQ, where the rate of work production is 
given by W = k{dE/dk). Then the NEQ work W done 
on the particle moving along a path q(T) for < t < t 
is given by 

This is rewritten in a matrix form as 



/3W[q] - i ^ drq^ ■ A ■ q 



(9) 



where 



-1/m 
k "f /m 



(3) 



A = /3H = 



Pk 




(10) 
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The system is assumed to be initially in EQ at /3 with 
k{0) = ki and will reach a final state with k{t) ^ kf, 
which is certainly far from EQ. In this situation, the JE 
states 



(11) 



where (• • ■) denotes the average over all possible paths 
q(T) and AJ^ is the free energy difference J-{kf) — J-{ki) 
at /3. The JE can be trivially derived from the Crooks 
relation Q 



Pp{W)=e" 



(12) 



where Pf{W) = {5{W - W[q]))F is the PDF for the 
work production W during the forward process with the 
change from ki to kf and vice versa Pr(W) for the reverse 
process. The JE and the Crooks relation can be proved 
for a general form of energy and external perturbation 
for the Langevin dynamics, if the system is initially in 
EQ at /3 [11|. However, the explicit expression for P{W) 
is not generally known. There are not many stochastic 
models which can be solved analytically for the PDF of 
fluctuating quantities. The breathing harmonic potential 
with a time-dependent force constant is not only analyt- 
ically tractable, but can also serve an appropriate model 
for the potential well in an optical tweezer or trap. 



III. PATH INTEGRAL FORMALISM WITH 
TIME-DEPENDENT FORCE 

The Fokker-Planck equation for a multivariate sys- 
tem with a linear drift force, known as the high di- 
mensional Ornstein-Uhlenbeck process, is solvable, i.e., 
the time-dependent PDF P(q, t) can be obtained ana- 
lytically [U . NEQ properties for this process were 
investigated in detail from the view of the circulating 
NESS current [l^ and the violation of the fluctuation- 
dissipation relation [4^ . Recently we have revisited this 
system in the light of the fluctuation theorem in the 
case that the (non-conservative) drift force docs not vary 
with time. The path integral formalism developed in that 
study can be extended to the present problem with time- 
dependent drift represented by the force matrix F(t) in 
Eq. © with k = k{t). 

To describe the NEQ fluctuations, it is convenient to 
introduce a path integral during time period t as 

/(qi,A;l(T))= j dqoPe,(qo;fc(0)) j i?[q] 

Xe^/o drL(^A)-XPWH+jl drl^.^ ^^^^^ 

The initial PDF for qo is chosen to follow the EQ 
Boltzmann distribution Peg(qo;fc(0)) as in Eq. 0, and 
/ £'[q](- • ■) denotes the integration over all possible paths 
connecting q(0) = qo and q(t) = qi for < r < The 
Lagrangian L is chosen to read 

i(q,q;F) = i(q+F.q)^-D-l.(q+F-q) . (14) 



The source term (J dr 1^ • q) is introduced for a later 
use. Note that the exponent of the integrand is at most 
quadratic in q. Hence the path integration can be com- 
puted exactly by Gaussian integrations. 

The quantity / is useful in calculating physical quan- 
tities of interest. For example, the PDF P(q, t) is given 
by 



P(q,t)=/(q,A-0;l(T) = O) 



(15) 



The PDF for the NEQ work production can be also 
calculated from /. First, we define a dimcnsionless quan- 
tity for the work as w = I3W for simplicity and introduce 
its generating function 

g{\) = (e-^'^^) = j dwe-^-^Piw) , (16) 

which can be obtained as 

(17) 



g{\)^ j dq/(q,A;l(r) =0) 



Note that the JE, Q{1) = exp[— /3AJ^], can be proven 
explicitly in this path integral formalism as well as 
the generalized Crooks relation as Qf{^)/Gb{^ ~ ^) ~ 
exp[— /3A7^] where F (R) denotes the forward (reverse) 
process. The PDF for the dimcnsionless work w is then 
obtained by the inverse Fourier transformation as 



(18) 



For an arbitrary functional .A[q(r)], one can also cal- 
culate its ensemble-averaged value from /. Defining the 
cumulant generating function as 



Z[l(r)] = Jdq /(q,A = 0;l(r)) 



one finds that 

(^[q]> = A 



61{t) 



mr)] 



(19) 



(20) 



1^0 



We will use this relation to calculate the cumulants of 
the work. 

The path integral, Eq. ([T^ . can be evaluated by using 
the methods developed in our recent study [l^]. Here, 
we will present the results without showing detailed cal- 
culation steps. 



A. Probability distribution function 

The PDF P(q, t) is given by 



P(q,t) ^ |det[27rA~^(t)] 



-1/2 



(21) 



where the kernel A(i) is a symmetric matrix, satisfying 
the differential equation as 



dA- 



dt 



2D - [F(t)A-i +A-iF^(t)] 



(22) 
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The formal solution is given by 

k-^{t)^2[ dT \J{t-t - T)m^ {t]t - t) 
Jo 

+ U(t;0)A-i(0)U^(t;0) (23) 

with the initial condition A(0) = /3H(0). Here the evolu- 
tion operator U is given by 



(24) 



TO 



where the subscript denotes the time-ordered product, 
and satisfies the differential equation 



d_ 

dt 



U{t;t') = -F{t)U{t;t') 



(25) 



with [J{t';t') = I (the identity matrix). 

In the absence of noises (D = 0), U{t; t') describes the 
deterministic evolution by q(i) — U(t;i') q(<')- When 
the force matrix is constant in time so that U(i,t') ~ 
g-it-t )F^ Qj^g (,g^j-^ ^YiQ integral in Eq. ([25)) and find 
the explicit solution for A(t) (see Eq. (21) in [13]). For 
a general time-dependent f{t), it is difficult to treat U 
analytically. However, Eqs. (^5)) and (pSj) can be solved 
very precisely by numerical integrations. 



B. Work distribution function 

The generating function for the work distribution in 
Eq. (jl7p involves the integration of the quantity / with 
nonzero A. The work W[q] coupled to A is also quadratic 
in q (see Eq. (|9])), hence the integration can be performed 
in the same way as was done for the probability distri- 
bution. After some algebra, one can derive 



Ina(A) 



dr Tr 



A-i(r;A)A(T) 



(26) 



where A — /3H in Eq. (|TOl) and A(t; A) is the modified 
kernel due to — A/3W in Eq. p^ . It is found to satisfy 
the nonlinear differential equation 



dA- 



dr 



2D - (FA-^ + A-^F^) - AA-^AA-i (27) 



where the initial condition is given by A(0; A) = /3H(0). 

This nonlinear differential equation can be solved eas- 
ily for A = and 1. The solution is A(r;0) = A(t) in 
Eq. while A(t; 1) = /3H(t). Interestingly, A(t; 1) 

corresponds to the kernel for P(q, r) in the quasi-static 
process. Inserting this into Eq. (|26p . we find 



which verifies the JE. 
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m 



= -/3AJ- , 
(28) 



It is not efficient to integrate Eq. (|27)) numerically, be- 
cause A sometimes becomes singular, as r increases in 
some range of A. Then, A~^ cannot be defined any more. 
Therefore it is more convenient to rewrite Eq. (j27p in 
terms of A(r; A) as 



dA 

d^ 



AA + (AF+ F^A) - 2 ADA 



(29) 



along with an equivalent and more effective expression 
for the generating function replacing Eq. (|26p as 



Ine(A) 



dr Tr 



F(T)-A(r;A)D 



1^ detA(t;A) 

2 '^detA(0;A) ' 
(30) 

A similar result has been found in the time-independent 
case [131 ■ Equations ((2^ and are ingredients for nu- 
merical study of the work production distribution P{w). 
Especially, the exponentially decaying tail behavior of 
P{w) is manifested by the divergence of ^(A), which 
turned out to be fully captured by the singularity in the 
logarithmic boundary term in Eq. pop . Thus we will 
focus on the behavior of det A{t; A) in the next section. 

One can observe that Eq. (|29)) becomes independent 
of /3 if A is scaled by (3. This proves that ^(A) as well 
as P{w) is independent of /?. Therefore, P{W) is simply 
equal to BP(w) with w = f3W. In the weak noise (large 
/?) limit [4lj, (4^, the tail behavior of P{w) for large |it; 
determines exactly and fully the work distribution P{W) 
except for a narrow central region, \W\ < /3"^. We will 
come back to this issue later. 



C. Cumulants of work production 

The cumulant generating function in Eq. (jl9p is found 



as 



^eM'^^I dr'F(r)-r(r,r')-l(r') 



(31) 



This form is expected because the Lagrangian is 
quadratic in q and the source field 1(t) is linearly coupled 
to q. The kernel r(T, r') is given as 

r U(r,T')A-i(r') ,r >t'; 

\ A-1(t)U^(t',t) ,t<t' ^ ^ 

and r(T,T') = r^(T',T). 

Using Eqs. ([201 and ([5T|) . one can express any func- 
tional of path q(T) in terms of \'{t, t'). For example, the 
first and the second cumulant of the work are given by 



H(T)r(r,T) 



(W) = ^j^^dT Tr 



(W')c-7t/ dr dr'Tr V[t,t')H{t')V^ {t,t')H{t 
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where (W^)^ = (W^) - (W)2. Note that = k and 
Hab = 0, otherwise. Then, the expressions become sim- 
pler: 



1 



(33) 



(W')c = ;^ / dT / dr' fc(T)MT')(rn(r,r')r-(34) 

One can also find higher-order cumulants in terms of 
V{t,t'), which are non-zero in all orders. This implies 
that the PDF P{w) should have a non-Gaussian form. 
The PDF shape will be discussed further in the next ses- 



D. Sudden change limit 

A sudden change is a rare case where one can calcu- 
late the work PDF exactly, even in the underdamped 
case. Suppose that the particle is in EQ under the har- 
monic potential with the force constant ki , and the force 
constant is changed abruptly to fc/ at time t = [4ll.l4^. 
If the particle is in a state q = {x,p)^ just before the 
change (t = 0"), its state still remains the same right af- 
ter the change {t = 0+) as well as the PDF P(q, t). Only 
change occurs in the potential energy, which results in the 
energy change AE = E{q, kf) - E{q, fcj) = ^{kf - k,)x^ 
for state q. Then the work production yV(q) ~ AE. 

As the initial distribution is given by the EQ Boltz- 
mann distribution, the PDF P{x) = J dp Peq{q;ki) = 
y^g^r/(2^g-/3fe.xV2^ Then, the PDF P{w) of the di- 
mensionless work w — j^W can be easily derived, using 
P{w)dw = 2P{x)\dx\, which yields that 



P{w) 



a > 



(35) 



with a = ki/{kf — ki) and 9{w) is the Hcavisidc step 
function. 

The generating function Q can be easily calculated, 
using Eq. as 



(36) 



which diverges at A = ki/{ki — kf) = —a as expected. 
The JE is also seen from ^(1) = (kf/ki)''^/'^. 

Our analytic formalism in the previous subsections can 
also reproduce G{X). The sudden change in the potential 
function can be studied by considering 

fc(r) = he{-T) + kfOir) , k{T) = [kf ~ h)5{T) . (37) 

Integrating Eq. (PH]) from t = 0^ to r = 0+, one get 
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A{Q+-\)^ A{Q--\) + \p{kf -k,)[ Q 



il~X)ki + Xkf 
1/m 



(38) 



where 1(0"; A) = /3i/(0-) is used. Then, from Eq. ([50]). 
one can easily reproduce the result in Eq. p6p . The 
cumulants of the work can be also easily calculated as 



{w) 



Iki 



and 



2k? 



(39) 



IV. ANALYSIS OF WORK DISTRIBUTION 

The analytic formalism developed in this paper is very 
useful to investigate the work distribution P{w) numeri- 
cally, in particular, its tail behavior, in contrast to direct 
numerical integration of the equations of motion where 
we always face with a statistics problem, becoming seri- 
ous in rare-event regions. In this session, we first present 
numerical data from the latter method to check fluctua- 
tion relations and get some insights on the nature of the 
work production distribution. Then, the tail behavior of 
P{w) is carefully examined by the former method. 

It is convenient to work with dimensionless variables 
in numerical calculations. Rescaling x = b^x, p = bpp, 

t = bti with 6, = y^, bp = bt ^ one 

finds that the dimensionless variables x and p satisfy the 
same equations of motion as in Eq. ([T]) with respect to 
the dimensionless time t with dimensionless parameters 
m — J ~ f3 ~ 1 and k = (^)fc. Hence, without loss of 
generality, we will set to = 7 = /3 = 1. The only indepen- 
dent parameter is the force constant k. We will consider 
a special case where fc is a time-independent constant, 
i.e., k{t) = ki{l + at) for convenience. 

First, we check the JE and the Crooks relation from 
direct numerical integration of the time-discretizcd equa- 
tions of motion. We adopt a notation A„ = X{t = i„) 
for a time-dependent quantity X{t) where i„ = nAt [n ~ 
0, 1, 2, • • •) are discretized times in unit of At. Then, the 
equations of motion are solved from the difference equa- 
tions 



{At)p„ 



i+1 = Pn - {At){pn + knXn) + \/2{Mj TJn 



where rjn are independent Gaussian-distributed random 
variables with zero mean and unit variance. An initial 
configuration qo = {xo,po) is drawn from the EQ dis- 
tribution of Eq. (O . The dimensionless NEQ work pro- 
duction Wn = PWn up to time i„ is evaluated from the 
recursion relation 

wn+i =wn + ^(a;'+i + xl)At (40) 

with wq = 0. Repeating the simulations Ng times, one 
can measure the PDF P{'w) and the generating function 
Cy(A) numerically. Note that the work production w is 
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FIG. 1. (Color online) (a) Pf{w) for the forward (F) process 
and Pr{w) for the reverse (R) process. Open symbols repre- 
sent e"'~'''2pB(-w). (b) ^^(A) for the F process and ^^(A) 
for the R process. Filled symbols represent the JE points. 
Open symbols represent e~'"'^C/^(l — A). Also shown with 
dashed lines are the generating functions obtained from the 
analytic formula in Eq. (|3Up . 



always positive for a > and negative for a < 0, inde- 
pendent of noise realizations. 

In simulations, we take At = lO^'^ and Ns — 10^. The 
force constant k{t) is taken to vary linearly from ki ~ 1 
to kf = A and fc^ = 4 to fc/ = 1, which will be referred 
to as a forward (F) and a reverse (R) process, respec- 
tively. Figure [TJa) shows the Pp (w) for the F process 
till i = 3 with a — 1 and Pr{w) for the R process with 
a = -1/4. We compare Pf{w) and e'"'^^-'^ Pr{~w) 

with (3 AT ~ 5 In ^ = In 2. They seem to overlap each 
other well (except for the region with very small P{w)), 
which supports the validity of the Crooks relation in 

Eq. dni). 

In order to examine the PDF in detail, we compute the 
generating function ^(A) = {e^^^'). These are plotted in 
Fig. mb). The JE states that ^(A = 1) = e'l^^^ where 
PAT = In 2 for the F process and — In 2 for the R process. 
Indeed, the numerical curves pass through the JE points. 
We also compare ^/_f(A) and e~^^-^QR{l — A), to check 
the generalized Crooks relation. For moderate values of 
A, both data align along a single curve. However, there is 
a slight discrepancy for large |A| where rare fluctuations 
with large values of |w| are important. This reflects a 
statistical uncertainty due to limited samplings. Even 
with Ns = 10^ samples, statistics is poor for those rare 
fluctuations. 

Now, we utilize the analytic results in Eqs. (|29)) and 
pop which arc free from statistical errors, in order to de- 
termine the tail part of P{w) precisely. In discretized 
times in unit of At = 10^'^, the nonlinear differential 
equation ((29)) for A(t; A) is solved with the initial condi- 
tion A(0;A) = /3H(0) and the integration in Eq. is 
performed numerically. We present the numerical results 
for the F and R processes (dashed lines) in Fig. [IJb). As 
expected, the previous simulation results deviate signifi- 
cantly from our new improved numerical results in rare- 




FIG. 2. Time evolution of det A(t; A) for the case with k(t) 
ki{l + at) with fci = 1 and a = 1 at various values of A. 



event regions. We checked that the relation Qf{^) = 
e~^^'^C/fl(l — A) is satisfied perfectly well with our new 
numerical results at all values of A. 

In fact, our numerical data in Fig. [ijb) shows that 
^(A) is divergent at a threshold Aq (F) and 1 — Aq (R) 
with Aq ~ —0.84713 < 0. The divergence occurs when 
detA(t;A) = as seen in Eq. (|30l). Figure [2] shows the 
time evolution of det A(t; A) at several values of A in the 
case with ki = I and a = I. To a given value of t, 
det A becomes smaller as A decreases and vanishes at a 
threshold Aq. One can solve the equation detA(t;A) = 
numerically to obtain the t-dependent threshold Aq. 
Figure [3] shows the numerical results for the system with 
ki ^ 1 and a ~ 0.5,1, and 2. The threshold depends 
on ki and a, and increases monotonically and converges 
to a finite limiting value Ag° ~ —1.39162, —0.81311, and 
-0.48110 in the t ^- oo limit. 

The singular behavior of Q{X) reveals the asymptotic 
behavior of the tail shape of P{w) for large \w\. Due 
to the generalized Crooks relation, it suffices to consider 
the F process with positive a (compression) . Figure [2] 
suggests that det A(t; A) is regular near A = Xo{t) so that 
one can write as det A ~ c(A — Ao(t)) with a positive 
constant c. Then, from Eq. ^(A) diverges as 

g{x)^{x-Xo{t))-'/' , (41) 

and its square-root singularity indicates that [l^] 

P{w) - i(;-l/2g-'|Ao(t)|.^ (42) 

for large positive if with the characteristic work wq = 
l/|Ao|. 

Our analytic formalism, combined with the numeri- 
cal analysis, predicts that there is the exponential tail 
in P{w) with the power-law prefactor with the expo- 
nent — 1/2. Note that the abrupt change of k (sudden 
change limit) also yields the same tail, as was shown in 
the preceding section. We test the tail shape by direct 
numerical integration of the equations of motion, using 
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t 



FIG. 3. (Color online) t-dependence of the threshold Ao for 
the F process with ki = 1 and a = 0.5, 1, 2. Multiple stepwise 
increases are observed in the insets showing the magnification 
of the curve with a = 1. 



lo' 




FIG. 4. (Color online) Rescaled PDF for the F process with 
ki = 1, kf = 4, and a — 1,2, 4, 8. The dashed line has a slope 
of -1/2. 

ki = 1, kf = 4, and various a — 1,2,4,8. For each 
case, the threshold Aq is obtained by solving the equation 
dct A(t; Ao) with fixed t = {kf - ki)/kia. In Fig. El 
the PDF P{w) multiplied with el^"!™ follows a power-law 
scaling for large w, which confirms the tail shape. Huge 
fiuctuations for large w are due to statistical errors in 
sampling rare events. 

The tail shape of P{w) in Eq. (|42|) is consistent with 
previous findings in the overdamped limit [4ll-l43|. It is 
also not surprising to find that |Ao(t)| decreases with t 
because one can easily expect that the work PDF should 
be more distributed (flatter) as t increases. However, the 
monotonic behavior of |Ao(i)l is not trivially smooth, but 
has an interesting repeating structure. 

In Fig. [3l we observe a stepwise change of |Ao(t)| in 
time, composed of a rather fast linear change followed 
by quite slow plateau-type change, which repeats itself 
but with decreasing size both in magnitude and time 



period, and finally converges to the limiting value of 
Ag°|. This implies that the exponential tail of P{'w) re- 
laxes into the limiting distribution via multiple (possi- 
bly infinitely many) fast-slow-type relaxation dynamics. 
These repeated fast-slow-type dynamics resemble multi- 
ple locking-unlocking dynamic transitions found in two- 
dimensional linear diffusion systems in the overdamped 
limit §4j. However, our case shows rather smooth 
crossovers between fast and slow dynamics, in contrast 
to sharp transitions with completely flat plateaus in 
Ao(i)| 13J]- '^^11 stepwise changes in our case 

as pseudo locking-unlocking transitions. In the mathe- 
matical language, we can not find any det A{t] A) curve 
tangential to the t axis (det A = 0) in Fig. [5J which pro- 
hibits a completely flat plateau, so no sharp transition is 
realized. 

It is easy to recognize that the oscillatory feature of 
det A in Fig. [5] evokes the stepwise change of |Ao(i)|. 
First, note that all det A curves show oscillatory wig- 
gles almost simultaneously in time and the oscillation 
frequency grows as t increases. So we can define a set of 
characteristic times {tf,t2, ■ ■ •) where all curves show a 
local minimum (-I-) or maximum (— ) simultaneously, at 
least, approximately. The oscillatory behavior is related 
to the increasing frequency of the harmonic oscillator 
caused by the increasing force constant k{t) = ki{l + at). 

Due to this oscillatory feature of the det A curves, one 
can easily figure out that the curves cross the t axis 
sparsely right after t = until t — , and densely 
for < t < , and so on. Therefore, Aq increases very 
fast during < t < and very slow during < t < , 
and this fast-slow relaxation dynamics repeats itself with 
increasing frequency. 

The underlying mechanism of these pseudo locking- 
unlocking transitions should be similar to one for the 
sharp transitions found in the two-dimensional linear dif- 
fusion systems [s^ . Only differences are the nature of the 
rotational current, which exists here only in the phase 
space of {x,p) and the time-dependent external force, 
which acts a role of the rotational driving force as well 
as the (time-dependent) anisotropic harmonic potential 
in the phase space. However, we could not find a sharp 
dynamic transition in our model with an arbitrary choice 
of parameters {ki, a). Recalling what we learned in [3^ . 
we guess that the anisotropy may be always small in our 
model, compared to the driving force magnitude, in or- 
der to avoid a sharp transition. For full understanding, 
however, a further investigation is necessary. 

For the overdamped one-dimensional case, we cannot 
have any rotational current, so the oscillatory behavior 
is expected to be absent at all, which is confirmed rigor- 
ously in the next session. Therefore, it can be concluded 
that the pseudo locking-unlocking transitions found in 
the underdamped case originate from the existence of 
the rotational current in the phase space. 
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V. OVERDAMPED LIMIT 

In the overdamped limit, the usual Fokker-Planck 
equation, replacing the Kramers equation, reads 



dP{x,t) _ d_ 
dt dx 



dx 



^[j-^k{t)x + {j(3)-'i-]P{x,t) . (43) 



Then one can use the analytic formalism developed so far 
for the Brownian dynamics by replacing D with (7/3) ^, 
F with 7~^fc, H with fc, and H with k. Then, the work 
generating function t/(A) is given by 



Ing(A)^ r..^^(^) 
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7/3 



1 A{t;X) 
— t: hi — , 

2 A(0;A) 



(44) 

where the scalar quantity A{t; A) satisfies a nonlinear dif- 
ferential equation 



^^ = (/?fc)A + 2^i-Ai2 
ar 7 7/3 



(45) 



with the initial condition A(0;A) I3k{0) = jSk,. We 
will set J — (3 ~ 1 without loss of generality. The 
overdamplcd limit is investigated in detail with different 
choices of /c(t). 



A. fc(r) = fci(l + ar) 

All the relevant informations are obtained from 
Eqs. (|44| and (|45|) . Unfortunately, the closed- form so- 
lution for A{t; X) or G{X) is not available. However, the 
highly accurate numerical solution is possible, which is 
shown in Fig. [SJa) for A{t;X) with ki = 1 and a = 1. 
As in the Brownian dynamics, it becomes zero at a t- 
depcndent threshold Aq. The threshold is plotted in 
Fig. [5]Jb). ^From the similar analysis for Eq. (|4ip . the 
PDF P{w) can be found to have the same tail shape as 
in Eq. (|42)) . Note that the particle dynamics does not 
display any oscillatory motion in the overdamped limit, 
hence A docs not either as seen in FiglSJa). Thus, the 
threshold Aq in Fig.[5]|^b) varies in time smoothly, showing 
no stepwise change at all. 

Engel and Nickelscn studied the same harmonic poten- 
tial problem in the overdamped limit (4ll . l42j . In their 
studies, they evaluated the path integral using the saddle 
point method in the low noise (large /3) limit. However, 
as pointed out in section Hill P{w) is independent of /3, 
so P{W) with w = PW can be exactly determined from 
the tail behavior of P[w) except for \W\ < /3~^. 

In order to calculate the cumulants of the work pro- 
duction, we need to evaluate the PDF kernel A{t), first, 
using Eqs. p3)) and (p4)) . As we do not need to deal 
with the time ordered product, we simply gets U (t, t') = 
J^,dTk{T)]. Then, one can explicitly calculate the 
cumulants for the work production. Using Eqs. (|23p and 





FIG. 5. (Color online) (a) Time evolution of A{t; A) in the 
overdamped case with fc; = 1 and a = 1 at various values of 
A. (b) t-dependence of the threshold Ao for the processes with 
fci = 1 and a = 0.5, 1. 2. 



(132]), we find 

A-\t) = 2 rciTe-'/>"'"^"'^+yfcrie-'/>"'("^ , (46) 
Jq 

and 

r{t,t')^e-I^'''^''^^^A-\t') 



where we use A ^(0) = k^ ^. 
For ^(t) = ki{l + ar), we find 



A~^{t) = k-'e 



1 „-2kit-kiat'' 



and 



where 



2kiT-\-kia.T 



g{t) = l + 2h dr e 
Jq 

The first and second cumulants are found as 



a 



2 Jo 



5(t) 



v2 /-t 



l\2 



2 Jo ^^Uo '^'^'^^'^''^'^^^'^'^ 
+ f dT'h{ry)g{rf\ , 



where 



ft-(r, T ) = e 



l\ _ ^-2fc,(T+r')-fe,Q(T2+T'^) 



(47) 



(48) 



(49) 



(50) 



(51) 



(52) 



(53) 



There are two extreme cases; quasi-static and sudden 
processes. For the quasi-static process, one can take the 
limit: a — >■ 0, t — >■ oo with a finite value of at ~ {kf — 
ki)/ki. Changing variable as u = kiaT/{kf — ki), one 
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can get an approximate result for the time integral. The 
important ingredient for the integration is 




cb^ 



du 



where c = (fc/ — I (kio) — > oo. Then, one can get 



kf\ a k'j — kf 



8 hkj 



(54) 



which agree with the results by Speck [43|. It indicates 
that the work distribution function is perfectly Gaussian 
centered around w = 13 A up to ©(oV and the non- 
Gaussianity starts to appear in 0{a^) |5l| . In the quasi- 
static limit (a — ^ oo), the work distribution function be- 
comes a delta function, as expected for EQ processes. 

For a sudden process, one can take the opposite limit; 
a — > oo, t — > with a finite value of = {kf — 
ki)/ki. Also using the same variable u, the integrand 

of du e'^" can be expanded in orders of c in the c — > 
limit. As a result, one can get 



{w) 



kf fci 



2h 

{kf - 
2fc? 



1 - ^^4^-^ U 0(a-^) 



3kia 
Ijkf-h) 
3q; 



+ 0{a-'^) (55) 



Note that {w) and (u'^)c are finite even for an instan- 
taneous change (a = oo), which agrees with the sud- 
den change limit for the underdamped case in Eq. (p9| . 
It is different from the case for the rigid wall moving 
with speed v in the v ^ oo limit, where we expect 
(w) 0. The difference is due to the distinctive 
situations. For the former, the collision occurs every- 
where with the harmonic potential, while for the latter 
the collision occurs only at the descending wall. The 
similarity lies in the non-trivial fluctuation around the 
average value. 



B. 



fc(r) = ki/{l + ar) 



With this specific form, one can find the closed-form 
solution for the work generating function ^(A). For a > 
0, the harmonic potential becomes flatter with time r > 
and the work w done on the particle is always negative. 
On the other hand, for a < 0, the harmonic potential 
becomes stiffer with time r (0 < r < l/|a|) and w is 
always positive. 

It is convenient to change variables as 



fx{u) = {l + aT)A{T;\) 
u=— ln(l -|- ar) . 



a 



(56) 
(57) 



(a) 



0.6 



0, 



(c) 



'okac-k) 
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FIG. 6. Curves representing the relation between ut = 
^ ln(l -I- at) and Aq = afciA()/2 for c < in (a), < c < fc^ in 
(b), and ki < c in (c). The values of (c, ki) axe taken to be 
(-1, 1), (1/2, 1), and (2, 1), respectively. 



The time-like variable u is monotonically increasing with 
T, starting from to oo for any nonzero a. ^^From 
Eq. ps)) . one obtains a differential equation for fx{u): 



dh 
du 



-2 [{fx 



(58) 



with K, = vA — c^, c = {2ki + q:)/4, and A = aki\/2. 
Note that n may be either positive real or pure imaginary 
depending on the range of A. In either case, the solution 
is given by 



fx{u) 



ki cos{2ku) + {cki — A) 
cos(2/s:ii) + [ki — c) 



sin(2Kii) 



sin(2K'ii) 



(59) 



with cos(ia;) ~ coshx and sin(ix) = isinha; for any x. 

With the solution for f\{u) or equivalently for A{t\ A), 
one can obtain the work generating function using 
Eq. It is useful to note that fx{u) c + 

ln[cos(2Ku) -I- {ki — c)sin(2Ku)/K]. After a straight- 
forward algebra, we find that 



a(A) 



cos(2KUt) 



cfcj-A sin(2KMf ) 
ki K 



(60) 



with Ut = i ln(l + at). 

The work PDF P{w) can be obtained by the inverse 
Fourier transformation of t/(A) in Eq. (|18p . Note that 
the generating function has an inverse square-root sin- 
gularity at a particular value of A = Ao(m) at which the 
denominator in Eq. (j60p vanishes. In fact, the singular- 
ity occurs when A{t;X) = or equivalently fx{u) = 0, 
as seen in Eq. (|44|) . The inverse square-root singular- 
ity in G{X) at A = Ap implies an exponential tail with a 
power-law prefactor [27[ in P{w) 



P{w) 



11/2 



j,Ao(«)u 



(61) 



in the uj — > — oo limit for a > (Aq > 0) or in the ui 
limit for a < (Aq < 0). 
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FIG. 7. versus a/ki. The dashed hnes are a flat straight 
hne of Ag° = 1 starting from the a/ki = 2 point (open circle), 
and of Ao° = starting from the a/ki = —2 point (filled 
circle) . 



The singular point satisfies the relation 



Ut 



tan 



Aq - cki 



(62) 



with Ao = akiXo/2. It should be understood that 



tan 



lix 



i tanh (x) and that < tan x < 



for a real x. Figure IHl shows the plots for the solution 
of Eq. ((62)) . where the divergence of ut is observed as 
Ao approaches the limiting value from above. Interest- 
ingly, the time dependence and the limiting value are 
very different depending on whether c < (a < — 2fci), 
< c < fcj {-2ki < a < 2ki), or k^ > c {a > 2k,). 
Especially, the limiting value Ag° = lim„^oo Ao is given 
by 



(a < -2h) 
, {-2k, <a< 2ki) 
1 , (a > 2h) 





{2ki+a 



(63) 



We present the plot of the Ag° as a function of a/ki in 
Fig.IZl 

There is an interesting symmetry of A§°(— a) = 1 — 
Ag°(a). This comes from the Crooks relation. The 
reverse protocol with respect to the forward protocol 
A:(r) = ki/{l + ar) should be given as fcr(r) = k{t — T) = 
kf/{\ + Utt) with Ur = ~akf/ki and kf = ki/{\ + 
at). If the system starts with the EQ distribution with 
kj.{Q) = kf, all results derived here can be also applied 
for the reverse process by replacing ki by kf and a by 
-akf/ki. Then, Eq. §^ gives us X^rM = Ag°(-a). 
The Crooks relation of Eq. HH) yields Aff = 1 - X^^ 
in the large w limit, which leads to our symmetry of 
X^{-a) = l-X^{a). 

We add a few remarks on the interesting a dependence 
of Ag°: (i) For a > 2ki, the tail shape of P{w) does not 
change with a with fixed Ajf = 1. When a is large 



enough, the harmonic potential flattens very fast. Then, 
the particle dynamics starting from the EQ distribution 
with ki would be rather localized and not fully relaxed 
into a flattened harmonic potential. So the fluctuation 
in w may be dominated by an initial transient behavior 
even in the long-time limit (t — > oo), independent of the 
detailed shape of k{t). The sudden change limit discussed 
in subsection IIII Dl corresponds to the a = oo limit with 
kf = 0, where Ag° = |fej/(fc/ - fc^)| = 1 from Eq. §5^ is 
consistent with the result for a > 2ki. Nevertheless, it is 
still quite remarkable to find Ag° = 1 for large but finite 
a. Similar features of initial-distribution dominance in 
the large deviation function in the long-time limit have 
been found in various different situations f30l - [33j . 

(ii) For I a I < 2ki, \Xoo\ decreases monotonically as |a| 
increases. This behavior is compatible with the common 
wisdom that the fluctuation gets stronger (longer tail in 
P{uj)) as the rate of the change in driving increases. 

(iii) When a < ~2ki (or c < 0), we obtain that Ajf = 
0. This implies that P(w) has a pure power-law tail in 
the positive- i(j region in the u— 7'0o(t— s-l/jal) limit. In 
this case, the driving is strong enough to generate huge 
fluctuations. 

Transient behavior of Ao is also investigated in two 
limits. In the short-time limit {ut 0), P(w) is expected 
to exhibit a delta function distribution centered at w = 0. 
This is confirmed by the solution of Eq. (|62|) given by 
Ao ~ ki/(2ut) or Ao — l/{aut) in the ut ^ limit. In 
the opposite limit {ut — > oo), Eq. ((62l) yields 



Ao 



4c' kj 4cMt 
ki-2c 



16^7 



(c<0) 
(c = 0) 
(0 < c < k,) 
(c = h) 



^ k,{2c-ki) + 



2c-fei 



-4(c-fci)ut 



(c > ki) 
(64) 

Note that the asymptotic behavior near A = A^ is very 
different, depending on the region. In terms of Aq and 
it is interesting to see a nontrivial power-law relaxation 
for a > 2ki {ki < c) such that Ao — 1 + ^^(1 + at)^^ with 
z = l- {2ki/a). 

The generating function also produces the cumulants 
of the work production by {w'^)c = dV^^^Q /d{—X)'^\^^Q. 
We focus on the mean value of the work, which is given 
by 



a 



fc.., + i(l-|)(l-e--"- 



(65) 



The quasi-static process corresponds to the limiting 
case where a — > 0, t — )■ oo with fixed at = {ki — kf)/kf. 
In this limit, we find 



{w)^-Un{l + at)^Un(^^ 



(66) 



11 



which agrees with Eq. ((54)) . For a sudden process, we 
take the opposite hmit where a — > oo, i — >• with the 
same fixed value of at in the above. Eq. (|65|) approaches 
(w) = {kf — ki)/2ki, which agrees with Eq. ([5^. 



VI. DISCUSSION 

The Brownian dynamics with mixture of position and 
momentum variables, having the even and odd parity re- 
spectively in time reversal, has not been fully scrutinized. 
In many literatures, the overdamped limit was taken for 
simplicity or as regards a light mass in experiment; other- 
wise the Kramers equation should be investigated, which 
is a nontrivial task. Putting the position and momentum 
on the same footing and introducing the singular diffu- 
sion matrix as in Eq. ([6]), we converted it into the usual 
Fokker-Planck equation, where the strict constraint on 
the position and momentum as S{x —p/m) is relaxed, so 
it is easier to handle. It is a well-known fact from the 
classical text books, but there are not many examples 
exploiting this method. In our study, we have proven 
this approach to be successful in finding the results ana- 
lytically and numerically beyond the overdamped limit. 

For the motion under the harmonic potential with a 
time-dependent force constant k(t), we have succeeded 
in examining the PDF of the work production rigor- 
ously, the most important quantity in NEQ fiuctuations. 
As a result, we have found the exponential tail with a 
power-law prcfactor in the PDF P{w) and |Ao(t)|, the 
characteristic constant in the exponential tail, to de- 
crease with time t showing a fine structure of infinite 
but not sharply-edged staircase. By comparing the mul- 
tiple locking-unlocking transitions (sharply-edged stair- 



case) found in the two-dimensional linear diffusion sys- 
tem (3JI , we call these rather smooth staircase as a man- 



ifestation of multiple pseudo locking-unlocking transi- 
tions. These pseudo transitions completely go away in 
the overdamped limit where no rotational current exists 
even in the phase space. We also consider some exactly 
solvable models in the overdamped limit and found an 
interesting power-law (not exponential) tail in P{w) for 
the case of rather fast compression (a < —2ki) with the 
protocol k{t) = ki/{l + at), which implies huge NEQ 
fiuctuations. 

The potential well in optical tweezers or an optical trap 
experiment is controlled by an external field, so can have 
the shape change in time due to a time-varying external 
field. Therefore, our study can serve a theoretical ba- 
sis for such experiments, since the potential well may be 
approximated to be harmonic. The perturbation theory 
might be exploited to investigate an anharmonic effect. 
Our recent study of the multi-dimensional diffusion dy- 
namics for a linear drift force (27l . Is^ in the overdamped 
limit can also be realized in such experiments. We sug- 
gest many interesting experiments to be carried out in 
this direction. 
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